\documentclass[a4paper,10pt]{article}
\usepackage{cite}
\usepackage{graphicx}
\usepackage{caption}
\usepackage{subcaption}
\usepackage[T1]{fontenc}
\usepackage[slovene]{babel}
%\usepackage[utf8x]{inputenc}
\usepackage{hyperref}
\usepackage{amssymb}
\usepackage{color}
\usepackage{pbox}
\begin{document}
    
\title{A Quantitative Model of Genetic Switches and Stability Analysis}
\date{\today}

\maketitle

\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}


%%%%%%%%%%%%%
%INTRODUCTION
%%%%%%%%%%%%%
\section*{Introduction}
In order to be able to better predict the behaviuor of \emph{in vivo} application of the genetic switches, we set up a quantitative model of a Transcription-activator like effector (TAL) based bistable switch. We apply experimental measurements for individual components and measure their behaviour when applied into a single system. Two different network topologies were quantified and compared. The essential comparable characteristics of a swith are stabilility (reaching a defined steady state), responsiveness (reaching the predefined steady state with the introduction of the corresponding inducer) and robustness (expressing equivalent qualitative behaviour under a wide range of conditions). The first model topology employs a simple, negative feedback loop. This setup is then enhanced by adding positive feedback loop and thus increasing all essential characteristics.   


%%%%%%%%
%METHODS
%%%%%%%%
\section*{Modelling}

In this section, the key modelling approaches wil be described. For a complete model and the parameter values of a bistable switch and bistable switch with a positive feedback loop, see Appendix A and B.

\subsection*{Inducible system}
A variety of well established models in the literature were examined, which resulted in choosing ordinary differential Hill equations model, presented in \cite{Tigges2009}. The former work presents a tunable synthetic mammallian oscillator, whose model was readjusted to present a switch struture. The host choice of mammalian cells and pristinamycin-based inducible system suits our model well. \\
The inducible system in both switches is built upon two different inducer molecules. In both models, we employ a pristinamycin inducible protein (PIP), which binds to its corresponding site upstream of the constitutive promoter (active state) and thus represses the transcription. After the introduction of pristynamicin, PIP is released from its bindind site (inactive state) and the transcription proceeds unrepressed. The described dynamics are modelled directly in the rate equation for the coresponding promoter: \\
\beq \frac{dPip_a}{dt} = Pip_p  \cdot  (\frac{Pc}{A_{Pc}+Pc}) \label{eq:pip} \eeq
\beq r_3 = k_3 \cdot f_v \cdot (\alpha_{Pc} + (1-\alpha_{Pc}) \cdot ( 1- \frac{Pip_a^{n_{Pc}}}{K_{Pc} + Pip_a^{n_{Pc}} })), \label{eq:pc} \eeq
where the parameters are as follows:  $Pip_a$ active PIP, $Pip_p$ inactive PIP, $Pc$ pristinamycin concentration, $A_{Pc}$ binding affinity, $r3$ effective transcription rate,  $k_3$ basal transcription rate, $f_v$ concentration scaling factor, $\alpha_{Pc}$ fraction of transcription rate when repressed, $Pip_a$ concentration of active PIP, $n_{Pc}$ Hill coefficient and $K_{Pc}$ bindind affinity. As the concentration of active Pip grows, the $(1- \frac{Pip_a^{n_{Pc}}}{K_{Pc} + Pip_a^{n_{Pc}} })$ factor approaches zero and the transcription rate approaches its repressed value $k_3 \cdot f_v \cdot \alpha_{Pc}$.  \\
Conversely, the rapalog inducible system works as follows: with the absence of the inducer, the transcription rate is equal to the basal minimal promoter rate and Rapalog inducible protein (RIP) is unbound (inactive state). After the introduction of an inducer, RIP changes to its active state and activates the transcription when binding to its corresponding binding site:
\beq \frac{dRip_a}{dt} = Rip_p  \cdot  (\frac{Rg}{A_{Rg}+Rg}) \label{eq:rip} \eeq
\beq r_4 = k_4 \cdot f_v \cdot (\alpha_{Rg} + (1-\alpha_{Rg}) \cdot \frac{Rip_a^{n_{Rg}}}{K_{Rg} + Rip_a^{n_{Rg}} }) \label{eq:rg} \eeq
where the parameters are as follows: $Rip_a$ active RIP, $Rip_p$ inactive RIP, $Rg$ Rapalog concentration, $A_{Rg}$ binding affinity, $r4$ effective transcription rate,  $k_4$ basal transcription rate, $f_v$ concentration scaling factor, $\alpha_{Rg}$ fraction of transcription rate, $Rip_a$ concentration of active RIP, $n_{Rg}$ Hill coefficient and $K_{Rg}$ binding affinity. As the RIP concentration grows, the transcription rate converges towards its activated value. Activation of the promoters was calculated as reporter activity in cells transfected with both the activator and reporter plasmids divided by reporter activity in cells transfected with the reporter plasmid only. On the average, the fold activation of the minimal promoter was measured to be about 40-fold. In similar manner as rapalog, the rapamycin and HecAct inducible system was modelled in the positive feedback loop model.\\

\subsection*{Transcription factors}
The core components of the system are Transcription activator-like effectors (TAL), combined with KRAB or VP16 small molecules to achieve transcription repression or activation respectively. The dynamics were modelled using Hill equations for repression and activation \cite{Alon2007a}. For example, the effect of a repressor on transcription rate is modelled as a factor in a rate equation:
\beq r_1 \cdot  \frac{1}{(1+\frac{(Tal57Krab_p)^{n_{57}}}{(s_{57} \cdot t_{57})^{n_{57}}})},  \label{eq:rep} \eeq 
The parameters used are: $r_1$ unrepressed transcription rate, $Tal57Krab_p$ concentration of TAL57KRAB protein, $n_{57}$ Hill coefficient, $s_{57}$ scaling constant, $t_{57}$ threshold concentration (denoting concentration producing half the maximal repression).
Intuitively, the larger the Tal57Krab concentration, the smaller the fraction of the remaining transcription rate, which leads to a state with lesser steady concentration level.\\
The fraction of the remaining transcription rate was determined from measurements of individual TAL97KRAB and TAL57KRAB contructs. The cells were transfected with variuos concentrations of repressor and reporter plasmids as well as reporter plasmids only. The fraction was calculated by comparing the relative luciferase activity of repressed and unrepressed reporter plasmids. The acquired data were furter interpolated to get an analitical function for repressor/reporter ratios and can be seen in Figure 1a. The acquired function is of the form:
\beq q_{57}(x) = a \cdot e^{b \cdot r} + c \label{eq:reptransfer},\eeq
where $a,b,c$ are fitted parameters, $r$ the input ratio of repressor/reporter plasmids and $q_{57}$ the resulting fraction of remaining transcription rate. The former is further converted in a scaling constant, $s_{57}$ and directly employed in equation \ref{eq:rep}. For example, if the input ratio of TAL57KRAB producing plasmids and plasmids containing TAL57 binding sites is 5.0, the remaning transcription rate will be about $10\%$ of the unrepressed rate. \\
In an analogue manner, the activation factor was modelled as follows:
\beq a_{97} \cdot \frac{Tal97:VP16_p^{n_{97}}}{(t_{97} + Tal97:VP16_p)^{n_97}} \label{eq:act}\eeq
The parameters used are: $a_{97}$ fold activation, $Tal97Vp16_p$ concentration of TAL97VP16 protein, $n_{97}$ Hill coefficient, $s_{97}$ scaling constant, $t_{57}$ threshold concentration (denoting concentration producing half the maximal repression). \\
To get a quantitative value for fold activation, we again extracted the knowledge from individual measurements for TAL97VP16 and TAL57VP16 constructs (Figure 1b). The fold activation factor was calculated  by comparing the relative luciferase activity of inactivated and activated reporter plasmids. Again an analitical function was acquired in the form of: 
\beq a_{97}(r) = c \cdot r ^{\frac{1}{b}}, \label{eq:acttransfer} \eeq
where $b,c$ represent interpolated parameters, $a_{97}$ resulting activation rate and $r$ the input ratio of activator/reporter plasmids. For example, if the input ratio of TAL97VP16 producing plasmids and plasmids containing TAL97 binding sites is $0.02$, the remaning transcription rate will be about $1000$ times the unrepressed rate. In the case of zztivators, because of the wide range and orders of magnitude of activation factors, we employ a constant mean value of the function over the relevant interval.

\begin{figure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/tal_krab_eff_curve.png}
                \label{fig:tal_krab_eff}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/tal_vp16_eff_curve.png}
                \label{fig:tal_vp16_eff}
                \caption{}
        \end{subfigure}
        \caption{(a) Transcription factor curve for Tal97Krab repressor complex (blue) and Tal57Krab repressor complex (red). The x-axis represents the input ratio of repressor-including plasmids versus plasmids including the corresponding binding sites. The y-axis shows the rate of repression, i. e. the remaining percentage of the original transcription rate during repression. (b) Transcription factor curve for TAL57/97:VP16 activator curve. The x-axis represents the input ratio of activator-including plasmids versus plasmids including the corresponding binding sites. The y-axis shows the rate of activation, i.e. the factor of increased transcription rate during activation.}
%TODO: fold-induction?
\label{fig:eff_curves}
\end{figure}

\subsection*{Stability Analisys}

To evaluate the performance of the bistable switch, stability analyses are performed on the derived model by exploring the system phase (state) space \cite{Gardner2000,Zakharova2011}. The state of the system is defined by a combination of values of the NEPTUN and mCITRIN reporter concentrations. Analogous to electronic digital circuits, we define LOW and HIGH states, using an empirically determined threshold of 100 nM. This leads to four possible states of the system, summarized in Table~1. \\
\begin{table}[] \footnotesize
    \begin{tabular}{|l|l|l|}
        \hline
    NEPTUN concentration & mCITRIN concentration & System state                            \\ \hline
        
    LOW & LOW & \emph{Ambigous} \\
    HIGH & LOW & \emph{NEPTUN State} \\
    LOW & HIGH & \emph{mCITRIN State} \\
    HIGH & HIGH & \emph{Ambigous} \\
    \hline
    \end{tabular}
    \caption{Possible states of the bistable switch system.}
\end{table}
During the simulation of the selected model, the a state is discovered in any time point, where the concentration of the reporters do not change, i.e.:
\beq \frac{dNep}{dt} = \frac{dMct}{dt} = 0\eeq
The former are represented in the phase space as fixed points. A system is said to be bistable if all of its fixed points are valid (non-amboigous) states. \\
The stability analysis was performed for intervals of different time values. For each parameter set, a simulation was run and the system performance was classified as one of the five possible, color coded behaviours, summarized in Table~2. 
\begin{table}[] \footnotesize
    \begin{tabular}{|l|l|l|}
        \hline
Color & Behaviour & Comment                            \\ \hline
    \color{green} Green  & Bistable  & The system exhibits both states correctly  \\
    \color{blue}  Blue & Monostable (NEPTUN) & The system exhibits NEPTUN state only  \\
    \color{red}  Red &  Monostable (mCITRIN) & The system exhibits mCITRIN states only  \\
    \color{yellow}  Yellow &  Ambigous & \pbox{20cm}{ Both reporter concentrations are either below or above \\ threshold simultaneously at the steady states} \\
    \color{black} Black  & Not valid & A steady state is not found \\
    \hline
    \end{tabular}
    \caption{Disjunctive classes of switch behaviours on a give time interval. The rules to determine the classification are evaluated for fixed points of the system on a given time interval.}
\end{table}














%%%%%%%%
%RESULTS
%%%%%%%%
\section*{Results}
\subsection*{Bistable switch}

%%%%%%%%
%BSW:RUN
%%%%%%%%
\begin{figure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bsw/BistableSwitch_s3_276.png}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bsw/BistableSwitch_s4_276.png}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bsw/BistableSwitch_s1_276.png}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bsw/BistableSwitch_s2_276.png}
                \caption{}
        \end{subfigure}
        \caption{Simulation of a bistable switch for $t = 14500$ min (a) Concentrations of the reporters. (b) Phase space of the system with state thresholds. (c) Concentrations of the inducers and active inducible proteins. (d) Concentrations of TAL97KRAB and TAL57KRAB proteins.}
\end{figure}

%%%%%%%%%%%%%
%BSW:ANALYSIS
%%%%%%%%%%%%%
\begin{figure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bsw/20120830202031446000_beef_bsw_01_py_plot_blur.png}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bsw/20120831130722273000_beef_bsw_02_py_plot_blur.png}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bsw/20120831170209643000_beef_bsw_03_py_plot_blur.png}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bsw/20120830202049207000_beef_bsw_07_py_plot_blur.png}
                \caption{}
        \end{subfigure}
        \caption{Stability analysis of a bistable switch with color code, described in Table~2. (a) State convergenge in the absence of inducers. (b) Pristinamycin effects on the system. (c) Rapalog effect on the system. (d) Hill coefficient variability effect.}
\end{figure}

To evaluate the performance of the switch, we ran a simulation for a long enough time interval for the system to exhibit both steady states (Figure 1). The introduction of inducer Pristinamycin and Rapalog put the system in NEPTUN and mCITRIN states respectively (Figure 1b). Upon the introduction of the inducer (Figure 1c), a sudden peak in corresponding reporter concentration arises. After inducer removal, the reporter and TAL concentration (Figure 1d) fall down by some degree, the reason being the reduced activity of the corresponding inducible promoter. The resulting concentration remains large enough to repress the opposing state and for the system to remain in the resired state after the inducer signal. Notice the difference in repressor effectiveness (Figure 1d), where TAL97KRAB binding constructs exhibit notably larger repression than TAL57KRAB, resulting in the system being on the edge of a defined and ambigous states (Figure 1b). Peak concentration at steady states for reporters were of the order of $1\mu M$. Time to switch from on state to another was of the order of 1000 min. \\
Using defined behaviour classes, characterized in Table~2, stability analisys revealed more stable switch characteristics. Without the introduction if an inducer on a given time interval it is not determined by design in which state the system will end up. The answer is given on Figure~3a. Various input ratios of plasmids containing TAL97KRAB and TAL97 constructs were evaulated. As predicted, the prevalence of a give plasmid resulted in system ending up in a correspondin reporter state. Large and equal input masses for both plasmids produced amboigus states. To ensure the proper functionality of the switch, instant introduction of an inducer is required, since large difference in input masses can result in the system not being able to switch between states (data not shown). \\
Furthermore, switching from one state to the other was characterized in terms in the duration and the time difference in the introduction of the inducers. The parameters were chosen so that the system exhibits one switch in a time interval of 7200 min with the corresponding inducer introduced at $t=0$ and for 500 minutes duration. The NEPTUN $\rightarrow$ mCITRIN change (Figure 3b) works quite differently than mCITRIN $\rightarrow$ NEPTUN (Figure 3c). In particular, pristinamycin must be present for at least 467 min to change state effectively, while the rapalog duration of as little as 60 min proved to be enough. For the system to exhibit both states in the time interval, a similar difference interval in removal/introduction of ~2000 min was required. \\
Various sources examine cooperativity as the condition for bistability \cite{Chatterjee2008, Malphettes2006}. The interval $[1.0,5.0]$ of Hill coefficient for TAL97KRAB and TAL57KRAB proteins were examined (Figure 3d). The correct behaviour of was reported for coefficients greater than ~2.3 for TAL97KRAB and greater than ~2.7 for TAL57KRAB respectively.





\subsection*{Bistable switch wih a positive feedback loop}
The addition of a positive feedback loop to the circuit proved beneficial to overall performance of the system (Figure 4). Notable are the absence of fluctuations on the inducer removal and one order of magnitude larger stable state concentration of reporters (Figure 4a). Consequently, the repression of the currently inactive constructs is much more effective, which is proven by non-noticable activity of repressed promoters. Fluctuations remain in the concentrations od TAL-KRAB and TAL-VP16 constructs, which does not affect the system state (Figure 4d). \\
Stability analysis confirms the improvement gained from the positive feedback loop (Figure 5). The absence of inducers again results in one of the states being exhibited, this time with notably smaller ambiguity area (Figure 5a). Due to better repressor performance an constitutive promoter leakage, the NEPTUN state is favourable for similar enough input masses. Again, prevalence in mCITRIN/TAL57KRAB input constructs result in the corresponding state prevailing. \\
Surprisingly, characteristics of the inducible system changed considerably. Pristinamycin minimal input duration improved, requiring only 160 min to perform a sucessful switch, which is $34\ \%$ the initial duration (Figure 5b). Rapamycin minimal input time was 494 min (Figure 5c). The in-between intervals were $2200$ min for the NEPTUN $\rightarrow$ mCITRIN change $1500$ min for mCITRIN $\rightarrow$ NEPTUN change, which is an improvement to the first model. \\
As predicted, significant difference in input masses of the plasmids could see problems arising in the transition from one state to another. Indeed, the input mass od TAL97KRAB containing plasmids above 210 ng forced the system into a monostable state (Figure 5d), while for masses smaller than 120 ng the system was forced into the other monostable state by introducing masses of TAL57KRAB containing plasmids larger than 280 ng. The notable difference in masses overhead is due to better repression rates of TAL97KRAB. \\
Furthermore, the system proved stable under whole intervals $[1.0,10.0]$ for TAL cooperativity coefficients, confirming the theory of positive feedback being able to make cooperativity a non-necessary condition for cooperativity. 
%TODO:
%http://www.annualreviews.org/doi/pdf/10.1146/annurev-biophys-050511-102240

%%%%%%%%%
%BFSW:RUN
%%%%%%%%%
\begin{figure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bfsw/BistableFeedbackSwitch_s4_730.png}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bfsw/BistableFeedbackSwitch_s5_730.png}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bfsw/BistableFeedbackSwitch_s1_730.png}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bfsw/BistableFeedbackSwitch_s3_730.png}
                \caption{}
        \end{subfigure}
        \caption{Simulation of a bistable switch for $t = 20000$ min (a) Concentrations of the reporters. (b) Phase space of the system with state thresholds. (c) Concentrations of the inducers and active inducible proteins. (d) Concentrations of TAL97VP16 and TAL57VP16 proteins.}
\end{figure}

%%%%%%%%%%%%%%
%BFSW:ANALYSIS
%%%%%%%%%%%%%%
\begin{figure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bfsw/20120831155929402000_beef_bfsw_01a_py_plot_blur.png}
                \caption{}
        \end{subfigure}

        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bfsw/20120831155933982000_beef_bfsw_02a_py_plot_blur.png}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bfsw/20120901002509723000_beef_bfsw_03a_py_plot_blur.png}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bfsw/20120831204103214000_beef_bfsw_04a_py_plot_blur.png}
                \caption{}
        \end{subfigure}
        \begin{subfigure}[b]{0.5\textwidth}
                \centering
                  \includegraphics[width=0.95\textwidth]{./img/bfsw/20120831155923907000_beef_bfsw_05a_py_plot_blur.png}
                \caption{}
        \end{subfigure}
        \caption{Stability analysis of a bistable switch with a positive feedback loop. Used is the color code, described in Table~2. (a) State convergenge in the absence of inducers. (b) Pristinamycin effects on the system. (c) Rapamycin effect on the system. (d) Plasmid dosage effect on NEPTUN$\rightarrow$mCITRIN state transition. (e) Plasmid dosage effect on mCITRIN$\rightarrow$NEPTUN state transition.}
\end{figure}











%%%%%%%%%%%%%%%%%%%%%%%%%
%A: BISTABLE SWITCH MODEL
%%%%%%%%%%%%%%%%%%%%%%%%%


\newpage
\section*{Appendix A: Bistable Switch model}

A schematic construction representation is depicted in Figure~6. The observed chemical species in the system are listed in Table~3, while the parameter values are found in Table~4. Ordinary diferential equations for the bistable switch model are Eq. 10-31.

\begin{figure}
  \includegraphics[width=120mm]{./img/Teamslovenia_bistable_sw.png}
  \caption{A bistable switch construction scheme.}
\end{figure}

\begin{table}[] \footnotesize
    \begin{tabular}{|l|l|l|}
        \hline
    Symbol & Chemical species                            \\ \hline
    $Rip_m$             & Rapalog binding protein (mRNA) \\
    $Rip_p$              & Rapalog binding protein  \\
    $Rip_a$              & Rapalog - Rapalog binding protein complex \\
    $Pip_m$              & Pristinamicyn binding protein (mRNA) \\
    $Pip_p$               & Pristinamicyn binding protein  \\
    $Pip_a$               & Pristinamicyn - Pristinamicyn binding protein complex \\
	  $Pc$ & Pristinamycin inducer \\
	  $Rg$ &  Rapalog inducer  \\
	  $Tal97Krab_m$ & Tal97Krab repressor protein (mRNA) \\
	  $Tal97Krab_p$ & Tal97Krab repressor protein \\
	  $Tal57Krab_m$ & Tal57Krab repressor protein (mRNA) \\
	  $Tal57Krab_p$ & Tal57Krab repressor protein \\
	  $Nep_m $ & NEPTUN reporter protein (mRNA) \\
	  $Nep_p $ & NEPTUN reporter protein \\
	  $Mct_m $ & mCITRIN reporter protein (mRNA) \\
	  $Mct_p $ & mCITRIN reporter protein \\
        \hline
    \end{tabular}
    \caption{The observed chemical species.}
\end{table}


\begin{table}[ht!]\footnotesize
    
    \begin{tabular}{|l|l|l|l|l|}
        \hline
        Symbol & Parameter & Value & Units & Source/reference                       \\ \hline
	$G_1$ & Plasmid 1 dosage & 190.0 & ng & Literature data \cite{Tigges2009} \\
	$G_2$ & Plasmid 2 dosage & 180.0 & ng & Literature data \cite{Tigges2009}\\
	$G_3$ & Plasmid 3 dosage & 180.0 & ng & Literature data \cite{Tigges2009}\\
	$G_4$ & Plasmid 4 dosage & 180.0 & ng & Literature data \cite{Tigges2009}\\
	$G_5$ & Plasmid 5 dosage & 180.0 & ng & Literature data \cite{Tigges2009}\\
	$G_6$ & Plasmid 6 dosage & 180.0 & ng & Literature data \cite{Tigges2009}\\
	$k_1$ & Plasmid 1 basal transcritiption rate  & 30.0 & mRNA $min^{-1}$ & Literature data \cite{Tigges2009} \\
	$k_2$ & Plasmid 2 basal transcritiption rate  & 30.0 & mRNA $min^{-1}$ & Literature data \cite{Tigges2009} \\
	$k_3$ & Plasmid 3 basal transcritiption rate  & 30.0 & mRNA $min^{-1}$ & Literature data \cite{Tigges2009} \\
	$k_4$ & Plasmid 4 basal transcritiption rate  & 30.0 & mRNA $min^{-1}$ & Literature data \cite{Tigges2009} \\
	$k_5$ & Plasmid 5 basal transcritiption rate  & 30.0 & mRNA $min^{-1}$ & Literature data \cite{Tigges2009} \\
	$k_6$ & Plasmid 6 basal transcritiption rate  & 30.0 & mRNA $min^{-1}$ & Literature data \cite{Tigges2009} \\
\hline
  $s_{97}$ & Tal97Krab scaling constant$^{\textrm{**}}$ & - & - &  Calculated from experimental data\\
  $s_{57}$ & Tal57Krab scaling constant$^{\textrm{**}}$ & - & - & Calculated from experimental data \\
  $n_{97}$ & Tal97Krab Hill coefficient & 3.0 &  - & Literature data \cite{Alon2007a}\\
  $n_{57}$ & Tal57Krab Hill coefficient & 3.0 &  - & Literature data \cite{Alon2007a}\\
  $t_{97}$ & Tal97Krab threshold concentration & $\tau$ & $nM$ & Adjusted by user defined state threshold  \\
  $t_{57}$ & Tal57Krab threshold concentration & $\tau$ & $nM$ & Adjusted by user defined state threshold \\
  
$d_{m}$ & mRNA degradation rate  & 0.0173 & mRNA $min^{-1}$  & Literature data \cite{Tigges2009}\\
	$d_{p}$ &  protein degradation rate & 0.0058 & nM $min^{-1}$ & Literature data \cite{Tigges2009}\\
	$t_{p}$ & translation rate & 0.2 & nM $min^{-1}$ & Literature data \cite{Tigges2009}\\
\hline
	$A_{Pc}$ & Pristinamicyn-Pip binding affinity & 1.0 & - & Literature data \cite{Tigges2009}\\
	$A_{Rg}$ & Rapalog-Rip binding affinity & 1.0 & - & Indentical to $A_{Pc}$\\
	$K_{Pc}$ & Promoter-Pip binding affinity & 1.0 & - & Literature data \cite{Tigges2009}\\
	$K_{Rg}$ & Rapalog-Rip binding affinity & 1.0 & - & Identical to $K_{Pc}$\\
	$\alpha_{Pc}$ & Basal activity of Pip inducible promoter & 0.085 & nM$s^{-1}$ & Literature data \cite{Tigges2009}\\
	$\alpha_{Rg}$ & Basal activity of Rip inducible promoter & 0.025 & nM$s^{-1}$ & Calculated from experimental data\\
	$n_{Pc}$ & Pip Cooperativity coefficient & 2.0 & - & Literature data \cite{Tigges2009} \\
	$n_{Rg}$ & Rip Cooperativity coefficient & 2.0 & - & Identical to $n_{Rg}$ \\
	%$h_{Pc}$ & Pristinimicyn half-life & 7.0 & min \\
	%$h_{Rg}$ & Rapalog half-life & 7.0 & min \\
	%$d_{Pc}$ & Pristinamicyn degradation rate$^{\textrm{***}}$ & - & nM $min^{-1}$ \\
	%$d_{Rg}$ &  Rapalog degradation rate$^{\textrm{***}}$ &-& nM $min^{-1}$ \\
\hline
	$f_{v}$ & Concentration scaling factor & $9.25 \cdot 10^{4}$ & - & Literature data \cite{Tigges2009, Batard2001} \\
	$\tau$ & Active state concentration threshold & 100 & $nM$ & Adjustable user defined state threshold \\
\hline	

        \hline
    \end{tabular}
    \caption{Model parameters and corresponding initial values.}
\end{table}


\beq r_1 = k_1 \cdot f_v \eeq
\beq r_2 = k_1 \cdot f_v \eeq
\beq r_3 = k_3 \cdot f_v \cdot (\alpha_{Pc} + (1-\alpha_{Pc}) \cdot ( 1- \frac{Pip_a^{n_{Pc}}}{K_{Pc} + Pip_a^{n_{Pc}} })) \eeq
\beq r_4 = k_4 \cdot f_v \cdot (\alpha_{Rg} + (1-\alpha_{Rg}) \cdot \frac{Rip_a^{n_{Rg}}}{K_{Rg} + Rip_a^{n_{Rg}} }) \eeq
\beq r_5 = k_1 \cdot f_v \eeq
\beq r_6 = k_1 \cdot f_v \eeq


%\beq \frac{dPc}{dt}  = - Pc  \cdot  d_{Pc}\eeq
%\beq \frac{dRg}{dt}  = - Rg  \cdot  d_{Rg}\eeq

\beq \frac{dPip_m}{dt} = G_5 \cdot r_5 - d_m \cdot Pip \eeq
\beq \frac{dPip_p}{dt} = t_p  \cdot  Pip_m - d_p \cdot Pip_p \eeq
\beq \frac{dPip_a}{dt} = Pip_p  \cdot  (\frac{Pc}{A_{Pc}+Pc}) \eeq

\beq \frac{dRip_m}{dt} = G_6 \cdot r_6 - d_m \cdot Rip \eeq
\beq \frac{dRip_p}{dt} = t_p  \cdot  Rip_m - d_p \cdot Rip_p \eeq
\beq \frac{dRip_a}{dt} = Rip_p  \cdot  (\frac{Rg}{A_{Rg}+Rg}) \eeq

\beq s_{97} = \frac{1.0}{\frac{1}{q_{97}}-1} \eeq
\beq s_{57} = \frac{1.0}{\frac{1}{q_{57}}-1} \eeq



\beq \frac{dTal97Krab_m}{dt} = G_1 \cdot r_1 \cdot  \frac{1}{(1+\frac{(Tal57Krab_p)^{n_{57}}}{(s_{57} \cdot t_{57})^{n_{57}}})} + G_3 \cdot r_3 - d_m  \cdot  Tal97Krab_m\eeq
\beq \frac{dTal97Krab_p}{dt} = t_p  \cdot  Tal97Krab_m - d_p  \cdot  Tal97Krab_p \eeq

\beq \frac{dTal57Krab_m}{dt} = G_2 \cdot r_2 \cdot  \frac{1}{(1+\frac{(Tal97Krab_p)^{n_{97}}}{(s_{97} \cdot t_{97})^{n_{97}}})} + G_4 \cdot r_4 - d_m  \cdot  Tal57Krab_m\eeq
\beq \frac{dTal57Krab_p}{dt} = t_p  \cdot  Tal57Krab_m - d_p  \cdot  Tal57Krab_p \eeq

\beq \frac{dNep_m}{dt} = G1 \cdot r1 \cdot  \frac{1}{(1+\frac{(Tal57Krab_p)^{n_{57}}}{(s_{57} \cdot t_{57})^{n_{57}}})} - d_m  \cdot  Nep_m\eeq
\beq \frac{dNep_p}{dt} = t_p  \cdot  Nep_m - d_p  \cdot  Nep_p\eeq

\beq \frac{dMct_m}{dt} = G1 \cdot r1 \cdot  \frac{1}{(1+\frac{(Tal97Krab_p)^{n_{97}}}{(s_{97} \cdot t_{97})^{n_{97}}})} - d_m  \cdot  Mct_m\eeq
\beq \frac{dMct_p}{dt} = t_p  \cdot  Mct_m - d_p  \cdot  Mct_p\eeq











%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%B: BISTABLE SWITCH WITH POSITIVE FEEDBACK LOOP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newpage
\section*{Appendix B: Bistable Switch model with a positive feedback loop}

A schematic construction representation is depicted in Figure~7. The observed chemical species in the system are listed in Table~5, while the parameter values are found in Table~6. Ordinary diferential equations for the bistable switch model are Eq. 32-56.

\begin{figure}
  \includegraphics[width=120mm]{./img/Teamslovenia_bistable_feedback_sw.png}
  \caption{A bistable switch with positive feedback loop construction scheme.}
\end{figure}


\begin{table}[] \footnotesize
    \centering
    \begin{tabular}{|l|l|l|}
        \hline
    Symbol & Chemical species                            \\ \hline
    $Hact_m$             & Rapamicyn binding protein (mRNA) \\
    $Hact_p$              & Rapamicyn binding protein  \\
    $Hact_a$              & Rapamicyn - Rapamicyn binding protein complex \\
    $Pip_m$              & Pristinamicyn binding protein (mRNA) \\
    $Pip_p$               & Pristinamicyn binding protein  \\
    $Pip_a$               & Pristinamicyn - Pristinamicyn binding protein complex \\
	  $Pc$ & Pristinamycin inducer \\
	  $Rg$ &  Rapamicyn inducer  \\
	  $Tal97Krab_m$ & Tal97Krab repressor protein (mRNA) \\
	  $Tal97Krab_p$ & Tal97Krab repressor protein \\
	  $Tal57Krab_m$ & Tal57Krab repressor protein (mRNA) \\
	  $Tal57Krab_p$ & Tal57Krab repressor protein \\
	  $Tal97:Vp16_m$ & Tal97:Vp16 activator protein (mRNA) \\
	  $Tal97:Vp16_p$ & Tal97:Vp16 activator protein \\
	  $Tal57Vp16_m$ & Tal57Vp16 activator protein (mRNA) \\
	  $Tal57Vp16_p$ & Tal57Vp16 activator protein \\
	  $Nep_m $ & NEPTUN reporter protein (mRNA) \\
	  $Nep_p $ & NEPTUN reporter protein \\
	  $Mct_m $ & mCITRIN reporter protein (mRNA) \\
	  $Mct_p $ & mCITRIN reporter protein \\
        \hline
    \end{tabular}
    \caption{The observed chemical species.}
\end{table}

\begin{table}[ht!] \footnotesize
    \begin{tabular}{|l|l|l|l|l|}
        \hline
        Symbol & Parameter meaning & Value & Units & Source/reference                          \\ \hline
	$G_1$ & Plasmid 1 dosage & 200.0 & ng & Literature data \cite{Tigges2009} \\
	$G_2$ & Plasmid 2 dosage & 200.0 & ng & Literature data \cite{Tigges2009}\\
	$G_3$ & Plasmid 3 dosage & 200.0 & ng & Literature data \cite{Tigges2009}\\
	$G_4$ & Plasmid 4 dosage & 200.0 & ng & Literature data \cite{Tigges2009}\\
	$G_5$ & Plasmid 5 dosage & 200.0 & ng & Literature data \cite{Tigges2009}\\
	$k_1$ & Plasmid 1 basal transcritiption rate  & 0.03 & mRNA $min^{-1}$ & Calculated from experimental data\\
	$k_2$ & Plasmid 2 basal transcritiption rate  & 0.03 & mRNA $min^{-1}$ & Calculated from experimental data\\
	$k_3$ & Plasmid 3 basal transcritiption rate  & 30.0 & mRNA $min^{-1}$ & Literature data \cite{Tigges2009}\\
	$k_4$ & Plasmid 4 basal transcritiption rate  & 0.03 & mRNA $min^{-1}$ & Calculated from experimental data\\
	$k_5$ & Plasmid 5 basal transcritiption rate  & 30.0 & mRNA $min^{-1}$ & Literature data \cite{Tigges2009}\\

\hline
	$a_{97}$ & Tal97Vp16 activation rate & -  & - & Calculated from experimental data\\
	$a_{57}$ & Tal57Vp16 activation rate & - & - & Calculated from experimental data\\
  $s_{97}$ & Tal97Krab scaling constant & - & - & Calculated from experimental data\\
  $s_{57}$ & Tal57Krab scaling constant & - & - & Calculated from experimental data\\
  $n_{97}$ & Tal97Krab/Vp16 Hill coefficient & 3.0 &  - & Literature data \cite{Alon2007a}\\
  $n_{57}$ & Tal57Krab/Vp16 Hill coefficient & 3.0 &  - & Literature data \cite{Alon2007a}\\
  $t_{97}$ & Tal97Krab/Vp16 threshold concentration & $\tau$ & $nM$  & Adjusted by user defined state threshold \\
  $t_{57}$ & Tal57Krab/Vp16 threshold concentration & $\tau$ & $nM$  & Adjusted by user defined state threshold \\
  
$d_{m}$ & mRNA degradation rate  & 0.0173 & mRNA $min^{-1}$ & Literature data \cite{Tigges2009}\\
	$d_{p}$ &  protein degradation rate & 0.0058 & nM $min^{-1}$ & Literature data \cite{Tigges2009}\\
	$t_{p}$ & translation rate & 0.2 & nM $min^{-1}$ & Literature data \cite{Tigges2009}\\
\hline
	$A_{Pc}$ & Pristinamicyn-Pip binding affinity & 1.0 & - & Literature data \cite{Tigges2009}\\
	$A_{Rg}$ & Rapamicyn-Hact binding affinity & 1.0 & - & Indetical to $A_{Pc}$\\
	$K_{Pc}$ & Promoter-Pip binding affinity & 1.0 & - & Literature data \cite{Tigges2009}\\
	$K_{Rg}$ & Rapamicyn-Hact binding affinity & 1.0 & - Identical to $K_{Pc}$\\
	$\alpha_{Pc}$ & Basal activity of Pip inducible promoter & 0.085 & nM$s^{-1}$ & Literature data \cite{Tigges2009}\\
	$\alpha_{Rg}$ & Basal activity of Hact inducible promoter & 0.025 & nM$s^{-1}$ & Calculated from exeprimental data\\
	$n_{Pc}$ & Pip Cooperativity coefficient & 2.0 & - & Literature data \cite{Tigges2009}\\
	$n_{Rg}$ & Hact Cooperativity coefficient & 2.0 & - & Identical to $n_{Pc}$\\
%	$h_{Pc}$ & Pristinimicyn half-life & 7.0 & min \\
%	$h_{Rg}$ & Rapamicyn half-life & 7.0 & min \\
%	$d_{Pc}$ & Pristinamicyn degradation rate$^{\textrm{***}}$ & - & nM $min^{-1}$ \\
% $d_{Rg}$ &  Rapamicyn degradation rate$^{\textrm{***}}$ &-& nM $min^{-1}$ \\
\hline
	$f_{v}$ & Concentration scaling factor & $9.25 \cdot 10^{4}$ & - & Literature data \cite{Tigges2009, Batard2001}\\
	$\tau$ & Active state concentration threshold & 100 & $nM$ & Adjustable user defined threshold\\
\hline
    \end{tabular}
    \caption{Model parameters and corresponding initial values.}
\end{table}



\beq r_1 = k_1 \cdot f_v \eeq
\beq r_2 = k_1 \cdot f_v \eeq
\beq r_3 = k_3 \cdot f_v \cdot (\alpha_{Pc} + (1-\alpha_{Pc}) \cdot ( 1- \frac{Pip_a^{n_{Pc}}}{K_{Pc} + Pip_a^{n_{Pc}} })) \eeq
\beq r_4 = k_4 \cdot f_v \cdot (\alpha_{Rg} + (1-\alpha_{Rg}) \cdot \frac{Hact_a^{n_{Rg}}}{K_{Rg} + Hact_a^{n_{Rg}} }) \eeq
\beq r_5 = k_1 \cdot f_v \eeq


%\beq \frac{dPc}{dt}  = - Pc  \cdot  d_{Pc}\eeq
%\beq \frac{dRg}{dt}  = - Rg  \cdot  d_{Rg}\eeq

\beq \frac{dPip_m}{dt} = G_5 \cdot r_5 - d_m \cdot Pip \eeq
\beq \frac{dPip_p}{dt} = t_p  \cdot  Pip_m - d_p \cdot Pip_p \eeq
\beq \frac{dPip_a}{dt} = Pip_p  \cdot  (\frac{Pc}{A_{Pc}+Pc}) \eeq

\beq \frac{dHact_m}{dt} = G_6 \cdot r_6 - d_m \cdot Hact \eeq
\beq \frac{dHact_p}{dt} = t_p  \cdot  Hact_m - d_p \cdot Hact_p \eeq
\beq \frac{dHact_a}{dt} = Hact_p  \cdot  (\frac{Rg}{A_{Rg}+Rg}) \eeq

\beq s_{97} = \frac{1.0}{\frac{1}{q_{97}}-1} \eeq
\beq s_{57} = \frac{1.0}{\frac{1}{q_{57}}-1} \eeq


%TAL97:KRAB
\beq \frac{dTal97Krab_m}{dt} = G_1 \cdot r_1 \cdot  \frac{1}{(1+\frac{(Tal57Krab_p)^{n_{57}}}{(s_{57} \cdot t_{57})^{n_{57}}})} \cdot a_{57} \cdot \frac{Tal57Vp16_p^{h_{57}}}{(t_{57} + Tal57Vp16_p)^{h_57}} + G_3 \cdot r_3 - d_m  \cdot  Tal97Krab_m\eeq
\beq \frac{dTal97Krab_p}{dt} = t_p  \cdot  Tal97Krab_m - d_p  \cdot  Tal97Krab_p \eeq

%TAL57:VP16
\beq \frac{dTal7Vp16_m}{dt} = G_1 \cdot r_1 \cdot  \frac{1}{(1+\frac{(Tal57Krab_p)^{n_{57}}}{(s_{57} \cdot t_{57})^{n_{57}}})} \cdot a_{57} \cdot \frac{Tal57Vp16_p^{h_{57}}}{(t_{57} + Tal57Vp16_p)^{h_57}} + G_3 \cdot r_3 - d_m  \cdot  Tal57Vp16_m\eeq
\beq \frac{dTal57Vp16_p}{dt} = t_p  \cdot  Tal57Vp16_m - d_p  \cdot  Tal57Vp16_p \eeq

%TAL57:KRAB
\beq \frac{dTal57Krab_m}{dt} = G_2 \cdot r_2 \cdot  \frac{1}{(1+\frac{(Tal97Krab_p)^{n_{97}}}{(s_{97} \cdot t_{97})^{n_{97}}})} \cdot a_{97} \cdot \frac{Tal97:VP16_p^{h_{97}}}{(t_{97} + Tal97:VP16_p)^{h_97}} + G_4 \cdot r_4 - d_m  \cdot  Tal57Krab_m\eeq
\beq \frac{dTal57Krab_p}{dt} = t_p  \cdot  Tal57Krab_m - d_p  \cdot  Tal57Krab_p \eeq

%TAL97:VP16
\beq \frac{dTal97Vp16_m}{dt} = G_2 \cdot r_2 \cdot  \frac{1}{(1+\frac{(Tal97Krab_p)^{n_{97}}}{(s_{97} \cdot t_{97})^{n_{97}}})} \cdot a_{97} \cdot \frac{Tal97:VP16_p^{h_{97}}}{(t_{97} + Tal97:VP16_p)^{h_97}}  + G_4 \cdot r_4 - d_m  \cdot  Tal97Vp16_m\eeq
\beq \frac{dTal97Vp16_p}{dt} = t_p  \cdot  Tal57Krab_m - d_p  \cdot  Tal97Vp16_p \eeq



%NEPTUN
\beq \frac{dNep_m}{dt} = G1 \cdot r1 \cdot  \frac{1}{(1+\frac{(Tal57Krab_p)^{n_{57}}}{(s_{57} \cdot t_{57})^{n_{57}}})} \cdot a_{57} \cdot \frac{Tal57Vp16_p^{h_{57}}}{(t_{57} + Tal57Vp16_p)^{h_57}} - d_m  \cdot  Nep_m\eeq
\beq \frac{dNep_p}{dt} = t_p  \cdot  Nep_m - d_p  \cdot  Nep_p\eeq

%MCITRIN
\beq \frac{dMct_m}{dt} = G1 \cdot r1 \cdot  \frac{1}{(1+\frac{(Tal97Krab_p)^{n_{97}}}{(s_{97} \cdot t_{97})^{n_{97}}})} \cdot a_{97} \cdot \frac{Tal97:VP16_p^{h_{97}}}{(t_{97} + Tal97:VP16_p)^{h_97}} - d_m  \cdot  Mct_m\eeq
\beq \frac{dMct_p}{dt} = t_p  \cdot  Mct_m - d_p  \cdot  Mct_p\eeq

\bibliography{Simulations}{}
\bibliographystyle{plain}

\end{document}
